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We consider the Bak-Tang-Wiesenfeld sandpile model on a two-dimensional square lattice of 
lattice sizes up to L — 4096. A detailed analysis of the probability distribution of the size, area, 
duration and radius of the avalanches will be given. To increase the accuracy of the determination 
of the avalanche exponents we introduce a new method for analyzing the data which reduces the 
finite-size effects of the measurements. The exponents of the avalanche distributions differ slightly 
from previous measurements and estimates obtained from a renormalization group approach. 
PACS number: 05.40,+j 






I. INTRODUCTION 



II. MODEL 
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Bak, Tang and Wicscnfcld |[J introduced the concept 
of self-organized criticality (SOC) and realized it with 
the so-called 'sandpile model' (BTW model). The steady 
state dynamics of the system is characterized by the 
probability distributions for the occurrence of relaxation 
clusters of a certain size, area, duration, etc . In the 
critical steady state these probability distributions ex- 
hibit power-law behavior. Using the concept of 'Abelian 
Sandpile Models' pi it is possible to calculate the static 
properties of the model exactly e.g. the height probabil- 
ities, height correlations, number of steady state config- 
urations, etc [J2|-|)|. However, the dynamical properties 
of the model, i.e. the exponents of the probability distri- 
butions, are not known exactly. Numerical simulations 
yield different values of the exponents depending on the 
considered system size and the used method of analyzing 
the data (see for instance pHPl). Recently Pietronero 
et al. [0| introduced a renormalization scheme which al- 
lowed them to estimate the avalanche exponents. An 
improvement of this renormalization scheme was given 
by Ivashkevich jl2| who obtained comparable results. 

We investigate the original Bak, Tang and Wiesenfeld 
model on large lattice sizes (L < 4096) and measured the 
probability distributions. Since the numerical investiga- 
tions of the BTW model by Manna || it is known that 
the obtained values of the exponents are affected by the 
finite size of the system. These finite size effects have 
to be taken into account in order to get the 'real' expo- 
nents. This has been done by extrapolation (L — ► oo) 
from data obtained for different L || . We could improve 
this method and are now able to measure the exponents 
of the infinite system directly thus avoiding any extrap- 
olation. In this way the accuracy of the obtained ex- 
ponents is increased significantly. We also address the 
question whether the BTW model and the related sand- 
pile models of Zhang |l3| and Manna [1J] belong to the 
same universality class. Finally, we discuss the assump- 
tion that the avalanche propagation can be described as 
a random walk. 



We consider the two-dimensional BTW model on a 
square lattice of size L X L in which integer variables 
hij > represent local heights. One perturbes the sys- 
tem by adding particles at a randomly chosen site hij 
according to 
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1 , with random (i,j). 
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A site is called unstable if the corresponding height hij 
exceeds a critical value h c , i.e. if hij > h c . Without loss 
of generality, we take h c — 4 throughout this work. An 
unstable site relaxes, its value is decreased by four and 
the neighboring sites are increased by one unit, i.e. 
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where the update is done in parallel. We assume open 
boundary conditions with heights at the boundary fixed 
to zero. 

System sizes from L = 64 to L — 4096 are investi- 
gated. Starting with a lattice of randomly distributed 
heights h e {0,1,2,3} the system is perturbed accord- 
ing to Eq. (fil) and Dhar's 'burning algorithm' is applied 
in order to check if the system has reached the critical 
steady state 0]. Then we start the actual measurements. 
All measurements are averaged over at least 10 6 non-zero 
avalanches except of the case L = 4096 where only 5 x 10 5 
measurements have been performed. We studied four dif- 
ferent properties characterizing an avalanche. In the fol- 
lowing we use the same notation as Majumdar et al. M. 
The total number of toppling events is called the size s 
of an avalanche. The number of distinct toppled lattice 
sites is denoted by Sd- Because a particular lattice site 
may topple several times the number of toppling events 
exceeds the number of distinct toppled lattice sites, i.e. 
s > Sd- The duration t of an avalanche is equal to the 
number of update sweeps needed until all sites are stable 
again. 
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FIG. 1. The probability distribution P(s) for different sys- 
tem sizes. The curves for L < 4096 are shifted in the down- 
ward direction. 

The linear size of an avalanche r is measured via the 
radius of gyration of the avalanche cluster. In the critical 
steady state the corresponding probability distributions 
should obey power-law behavior characterized by expo- 
nents t s , Td, Tt and r r according to 

Ps(s) ~ s- T % (4) 

Pd(sd) ~ s d ~ Td , (5) 

Pt(t) ~ t- r \ (6) 

P r {r) ~ r-^. (7) 

III. SIMULATIONS AND RESULTS 

Fig. [| displays the obtained results for the distribu- 
tion P s (s) for different system sizes. A power-law fit 
to the straight portion of these curves yields the expo- 
nents t s (L). Fig. D shows a plot of the exponents t 8 (L) 
vs. 1/lnL. It is seen that for L > 128 the exponents 
obey the finite-size behavior 



,W = 



const 
lnL 



(8) 



as suggested already by Manna |6|] . The extrapolation to 
L — > oo yields the value of the exponent T Si00 = 1.247. 
The probability distributions Pd(sd), Pt(t) and P r (r) are 
analyzed in the same way with the result Td,oo = 1-258, 
7t,oo = 1-405 and tv i00 — 1.588, respectively. All ex- 
ponents are slightly larger than those obtained from ear- 
lier simulations by Manna who considered smaller system 
sizes and had less statistics 101 . 



FIG. 2. Determination of the exponent r using the extrap- 
olation (Eq. Eh. 



However, these values of the exponents are not very 
accurate. Namely, a crucial point in this analysis is the 
extension of the fit region in each distribution P s (s,L). 
Changing it slightly different exponents are obtained. 
This uncertainty in the determination of the exponents 
t(L) can be estimated to be at least of the order of ±0.01. 
Taking then the propagation of these errors into account 
we can estimate the uncertainty in the determination of 
the extrapolated value Too to be of the order of ±0.06 
which is mainly due to the large distance of the mea- 
sured values from the vertical axis (see Fig. g). Thus it 
is in principle not possible to obtain the exponents of the 
BTW model with high accuracy by a simple extrapola- 
tion of the exponents via Eq. m). 

However, it is possible to improve the determination 
of the exponents not by using Eq. (||) for an extrapo- 
lation but for a direct determination of t^. Consider 
for this purpose two probability distributions P(s, L\) 
and P(s, L2) corresponding to different system sizes with 
L\ > Li- If Eq. (ph describes the finite-size behavior of 
the exponents t s correctly the probability distribution 
(Eq. [|) for a given system size L behaves as 



P(s,L) 



(9) 



Thus, the exponent r St00 can be determined directly by a 
power-law fit of the function H(s, L\, L 2 ) which is defined 
as 
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In Fig. H _ff(s,Li,L 2 ) is plotted for various system 
sizes L\ and L 2 . A nice property of this function is that 
in contrast to the probability distribution the cut off of 
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FIG. 3. The function H(s,Li,Li2) for different pairs L\ 
and 1/2- The curves are shifted with increasing system sizes 
in the downward direction. The solid lines correspond to a 
power-law fit. The obtained values of the exponent t 3 are 
listed in Table §. 



the power-law behavior at large values of s is now very 
abrupt. We apply this analysis to all four distributions 
and the resulting exponents are listed in Table |. The val- 
ues of the exponents t Sj00 , T tj00 and T rj0 o (except of the 
case L 2 = 128, Li = 256) fluctuate around their mean 
values given by r S!OC = 1.293±0.009, r t)OC = 1.480±0.011 
and T r ,oo = 1.665 ± 0.013. Only the exponent Td )00 dis- 
plays a significant L-dependence. A possible origin of 
this L-dependence is that Eq. (g) does not describe cor- 
rectly the finite size behavior of Td and that one has to 
add corrections to it. However, the data suggest that 
this additional L-dependence vanishes for large system 
sizes and therefore the exponent saturates in the vicin- 
ity of t^oc ~ 1.33. Note that the mentioned error-bars 
describe only the statistical error. Due to the systematic 
errors the real error-bars are slightly larger. 



IV. DISCUSSION 

Despite of their different toppling rules it is supposed 
that the BTW model, Zhang's model Q and Manna's 
Two-State model ]14| belong to the same universality 
class, i.e. they should be characterized by the same ex- 
ponents. Pietronero and co-workers [O] addressed this 
question by a renormalization group approach and found 
that the BTW model and Manna's Two-State model be- 
long to the same universality class. Different results were 
obtained by Ben-Hur and Biham p0| who found different 
values for the two models. 

In Table |H| we compare our results with the expo- 
nents of the Zhang and the Two-State model obtained 
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from recent investigations on comparable lattice sites 
Jl5]| . Within the error-bars the BTW and the Zhang 
model displays the same exponents. The differences of 
the exponents Td and r r of the BTW and Manna's model 
can not be explained by the error-bars and thus we con- 
clude that both models do not belong to the same uni- 
versality class. But it is remarkable that both models 
display nearly the same duration exponent Tt and espe- 
cially that Tj w |. We assume that the value r t — | 
is a common feature of many sandpile models caused by 
an analogy of the avalanche propagation and a random 
walk, which we will discuss now. 

The number of critical sites n(t) at a given update 
(time) step t can be considered as a random walker. 
Starting with n(t = 0) = 1 the avalanche performs a ran- 
dom walk n(t — 0) — » n(t = 1) — > n(t = 2) — ► ... with the 
transition probabilities p(n,n'). The avalanche ceases to 
exist if the random walk returns to the origin (n — 0). In 
the simplest case the transition probabilities are homoge- 
neous p(n,n') = p(n — n'), symmetric p(An) = p(— An) 
and the random numbers An are uncorrelated. Then the 
avalanche probability distribution Pt(t) is given by the 
probability Pf irat roturn (t) that a random walker with initial 
value n(t = 0) = 1, with certain transition probabilities 
for increasing, decreasing and maintaining n, returns for 
the first time to its starting point in t steps, which scales 
as HI 



P fl rst Mto n(t)~r 



(11) 



Certain sandpile models are solved by an exactly map- 
ping of the avalanche propagation onto a simple random 
walk |l|,|| . 

In Fig7T3 we present the number of critical sites vs. 
update steps of a certain avalanche of the BTW model. 
The probability distribution p(An) and the correspond- 
ing correlation function 

TABLE II. Values of the exponents r 3 , Td, Tt and r r for the 
BTW model, Zhang's model and Manna's Two-State model. 
Due to a finite curvature of the probability distribution the 
duration exponent Tt of the Zhang model can not be deter- 
mined in the usual way JL5|. 



Model 


Ts 


Td 


Tt 


T r 


BTW 


1.293 


1.330 


1.480 


1.665 


Zhang 


1.282 


1.338 




1.682 


Manna 


1.275 


1.373 


1.493 


1.743 
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FIG. 4. The avalanche propagation as a random walk. The 
number of critical sites n(t) is plotted against the update 
(time) steps for a certain avalanche of duration t = 289. 
Starting from n(t — 0) = 1 the avalanche stops if the random 
walker returns to the origin for the first time. 
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( An(t) An(t + At) ) 
(An*] 



(12) 



are shown in Fig. 0. The probability distribution p(An) 
has to be symmetric in order to make sure that the ran- 
dom walk is recurrent, i.e. the probability that it ever re- 
turns to the origin is one jig]. The distribution displays 
asymmetries only for finite system sizes. A detailed anal- 
ysis (not shown) yields that the third central moment of 
the distribution p(An) tends to zero with diverging sys- 
tem size L indicating that pi,— >oo(An) is symmetric. 

The correlation function C(At) is sharply peaked at 
At = but there are small oscillations for small values 
of At. Therefore, the second requirement for Eq. (|ll]) to 
be valid, uncorrelated steps An, is only fulfilled approx- 
imately. This oscillating behavior is caused by the used 
parallel update process. Since toppling occurs at a given 
time step in one sublattice only the update algorithm 
switches in sequently time steps between the two sublat- 
tices. The alternating correlation function indicates that 
the correlations within one sublattice differs from the cor- 
relation between the two sublattices. Thus, compared to 
the exact solved sandpile models J17yi8| where the cor- 
relation functions are simply given by a <5-function the 
correlations of the BTW model are more complicated. 
But since these oscillations at small At have small am- 
plitudes we suggest that the avalanche propagation may 
be described as a random walk and that the exponent of 
the duration is r t = |. 

Scaling relations for the exponents T s ,Td, Tt and r r can 
be obtained if one assumes that the size, area, duration 
and radius scale as a power of each other, for instance 



FIG. 5. The probability distribution p(An) (see in- 
set) and the correlations C(At) for different system sizes 
(L = 64, 128,256). 



for the duration t of an avalanche and its radius r. The 
relation P t (t)dt = P r (r)dr for the corresponding distri- 
bution functions leads to the scaling relation 



llr 



T r — 1 
Tt-1' 



(14) 



The exponents jdr, Irs, Isd etc are defined in the same 
way. The exponent jt r is usually identified with the dy- 
namical exponent z and using a momentum-space analy- 
sis of the corresponding Langevin equations Diaz-Guilcra 
showed that the dynamical exponent of the BTW and 
Zhang's model is given by z = (d + 2)/3 |jlj. On the 
other hand one concludes from the compactness of the 
avalanche clusters that jdr = 2. Thus one gets two scal- 
ing relations for the exponents Td, r r and r t and using the 
result that Tt = § the exponents of the probability dis- 
tribution of the radius and the area are given by r r = | 
and Td = §• These values are in good agreement with our 
numerical results and we would suggest that they are the 
exact exponents of the BTW model. 

Majumdar and Dhar assumed that the size and the 
area of an avalanche fulfill the relation 



s d n c 



(15) 



where n c is the number of topplings at the site initiating 
the avalanche. If this equation holds the exponents r s 
and Td have to fulfill the relation t s = 2 — l/r^. Using 
Td = o from above we obtain r s = | which is well outside 
the error-bars of our numerical result r s = 1.293. Thus 



yltT 



(13) 



we conclude that the assumed relation Eq. (15) does not 
describe the real scaling behavior. 

Due to a lack of a scaling relation which connects t s 
with the other exponents the exact value of the expo- 
nent To is still unknown. Even a numerical determination 
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FIG. 6. The conditional probability distribution p(s\sd) 
The arrow marks the corresponding expectation value. 



of the exponent j s d yields useless results. The relation 
which defines j s d implies the assumption that the condi- 
tional probability distribution p(s\s ( j) is strongly peaked 
so that the expectation value E(s\sd) scales with the area 
Sd- Measurements of the conditional probabilities show 
that this is not the case for p(s\sd) (see Fig. g). The dis- 
tribution displays an asymmetric shape which violates 
the above assumptions. 

A similar analysis of Manna's Two-State model yielded 
that the dynamical exponent is given by z rs § resulting 
in TV = |, T d = % H 



The BTW model and the Two- 
State model belongs to different universality classes. 



V. CONCLUSIONS 

We studied numerically the dynamical properties of 
the BTW model on a two-dimensional square lattice and 
measured for large system sizes (L < 4096) the avalanche 
probability distributions. We introduced a new analy- 
sis to minimize the finite-size effects and determined the 
avalanche exponents with an improved accuracy. Our 
numerical results are consistent with the values r t = f , 
r r = | and Td = I which we consider to be the ex- 
act exponents of the BTW model. We discussed the 
possibility that these values are caused by an analogy 
of the avalanche propagation and a random walk pro- 
cess. Further work has to be done to check this assump- 
tion. Recently Ivashkevich |12] improved the renormal- 
ization group approach for sandpile models proposed by 
Pietronero et al. uM. Both calculation yields the expo- 
nent Td ~ 1-25 significantly smaller than our numerical 
estimates. 
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